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The motion of contaminant particles through complex environments such as fractured rocks or 
porous sediments is often characterized by anomalous diffusion: the spread of the transported 
qq . quantity is found to grow sublinearly in time due to the presence of obstacles which hinder particle 

■ migration. The asymptotic behavior of these systems is usually well described by fractional diffusion, 

which provides an elegant and unified framework for modeling anomalous transport. We show 
£SJ . that pre-asymptotic corrections to fractional diffusion might become relevant, depending on the 

microscopic dynamics of the particles. To incorporate these effects, we derive a modified transport 
C$ ' equation and validate its effectiveness by a Monte Carlo simulation. 
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I. INTRODUCTION 

A diffusion process is called anomalous when the variance of the transported quantity, (x 2 (t)}, grows in time as 
t a , with a larger (superdiffusion) or smaller (subdiffusion) than one. A general approach to the analysis of stochastic 
transport is based on the continuous time random walk (CTRW), a probabilistic model where the motion of a walker 
in a medium is interpreted as a series of jumps of random lengths, separated by random waiting times [H [H S 0] ■ The 
theory of CTRW with power-law probability distribution functions (pdf 's) for the waiting times was first introduced 
in physics in the late 1960s to describe anomalous diffusion of charge carriers in amorphous semiconductors [1, [1,0, Hi- 
lt was later applied with success to a broad spectrum of physical systems (see, for example, Refs. Q and Q for a 
C$ . review). 

As a paradigmatic example of anomalous diffusion we consider here the case of contaminant transport through 
heterogeneous materials, such as porous sediments or fractured rocks. In this system the microscopic dynamics 
of each particle is assumed to be hindered by obstacles, dead ends, trapping events due to the interaction with 
the surrounding environment or abrupt changes in the flow field. The macroscopic effect is that the spread of the 
migrating particles plume grows sublinearly in time, thus resulting in subdiffusion f9l [Tol [Til IT2I Il3l [hil [TH flfl \vi\ . 
These assumptions are substantially corroborated by many experimental observations [l(J, [TJ, [l2|, HH, EDI ■ 

In general, the CTRW equations do not allow for closed-form analytical solutions. However, if suitable first-order 
approximations are introduced, a fractional diffusion equation can be formally derived from CTRW and analytical 
solutions can be obtained: in this respect, the fractional diffusion equation represents an asymptotic subset of CTRW 

A complementary approach to the modeling of stochastic transport is provided by Monte Carlo simulation, a particle 
tracking method which allows the fate of each walker to be followed by sampling the jump lengths and the waiting 
times from a suitable distribution. Monte Carlo simulation has been widely adopted as a natural tool to obtain 
reliable and accurate estimates of the CTRW solutions [ijj] . 
' In this paper, we explore the range of validity of the asymptotic fractional diffusion equation. This topic has been 
previously examined e.g. in Refs. [23[ and [24J. With the help of Monte Carlo simulations, building on the discussion 
in Refs. [l7| and (25|, we show that pre-asymptotic corrections to the fractional diffusion equations play a significant 
role, depending on the microscopic dynamics of the particles [42j]. Neglect of these corrections would lead to gross 
errors in the estimation of model parameters from measured data. To overcome this problem, we derive a modified 
equation which incorporates these effects. 

The paper is organized as follows. In Sec. [IT] we review the CTRW model and the asymptotic fractional diffusion 
equation. In Sec. IIHI we show that for a given range of a critical parameter the fractional diffusion equation might 
not properly characterize the actual spread of the transported quantity as derived from CTRW (through a Monte 
Carlo estimate). We discuss how to modify the fractional diffusion model so to include the contribution of pre- 
asymptotic corrections and validate the proposed results by comparing them to Monte Carlo simulations. Conclusions 
are discussed in Sec. IIV1 
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II. FROM CTRW TO FRACTIONAL DIFFUSION 



We briefly review here the main assumptions of the CTRW model. Consider a walker whose stochastic trajectory 
in a medium is modeled as a series of random jumps separated by random waiting times, during which the walker 
stays in the previously reached position. The associated pdf P(x, t) represents the probability density of the walker 
being at x at time t and is usually called the propagator of the process. For contaminant transport, P(x, t) represents 
the normalized particle concentration. The properties of the propagator depend on the jump lengths pdf A(x) and the 
waiting times pdf w(t). Once these distributions are assigned, it is possible to derive a probability equation for P(x, t), 
the master equation of the CTRW model 0,0,0,11] It can be shown Q that the master equation in {x, t} space can be 
recast into a simple algebraic relation for the Fourier and Laplace transformed propagator P{k,u) = C {!F{P(x, t)}} 

m 

V ' ; u l-\{k)w{uY v ; 

with the walkers starting at x ~ at time t = 0. 

Within the CTRW scheme the long waiting times which characterize the trajectory of a contaminant particle moving 
through heterogeneous materials are taken into account by assuming that w(t) has a power-law tail w(t) ~ 
a < 1, when t — > +oo. Therefore, w(t) lacks a characteristic time scale (the average of the distribution is infinite) and 
anomalously long waiting times have a non-negligible probability of being sampled. In the following we will assume 
a Pareto pdf of the form 

w(t) = a^r 1 "", (t > t) (2) 

where r is a time constant. This specific form for w(t) has been chosen so that Monte Carlo sampling by inverse 
transform is straightforward: it can be shown that the short time behavior of the pdf is of minor importance and 
only the tail matters in determining the properties of w(t) 27]. The jump length pdf A(a;) is usually assumed to be 
a Gaussian distribution, so that the jumps have a typical scale (say the standard deviation of the distribution) and 
extreme events, that is, jumps whose length are much larger than the standard deviation, are very improbable. 

In general, no analytical solution is known for Eq. ([1]), because the required inverse transforms are usually not 
trivial. However, exact results can be obtained for P{x, t) sufficiently far from the origin, in the diffusion limit 
\x\ — » +oo, t — * +oo, that is, k — > 0, u — » (see, for example, Ref. [2(| for a detailed discussion). 

If w(u) = £{w(t)} is expanded in the long time limit ut <C 1 and truncated to the first non-constant term in 
u, we find w{u) ~ 1 — c a (uT) a , where c a = T(l — a). More generally, it can be shown that any pdf of the kind 
w(t) ~ t~ 1 ~ a would lead to the same exp ans ion w{u) in the transformed space; the specific value of the constant c a 
depends on the details of the function @, H3]- I n contrast, if X(x) is a Gaussian with mean \i = and variance 2a 2 , 
then A(fc) = J- {X(x)} becomes A(fc) ~ 1 — k 2 a 2 in the diffusion limit ka <C 1. Any jump pdf with finite variance 
would asymptotically lead to the same result Q. 

We now introduce the Riemann-Liouville fractional derivative operator (28l. l29l. |30| 9q7"' which is defined by its 
action on a sufficiently well behaved function /: 

The Riemann-Liouville operator is an integro-differential operator involving a convolution integral with a power-law 
kernel. In particular, the operator 8q" may be thought as a generalization of the usual multiple integral of integer 
order. From its definition, <9^~" has a very simple representation in £-space: £{9(J~f/(f)} = u~ a f(u), Va > 

SiliiliS. 

If we substitute w(u) and A(fc) in Eq. (Q} and make use of the definition ([3]), we can formally write the propagator 
([1]) in the {x, f] space: 

8 f) 2 



-P^t) = D a d^ a —P{x,t), (4) 



where D a = a 2 /c a r a may be thought as a generalized diffusion coefficient Qj. Equation ((4]), which is called the 
fractional diffusion equation due to its fractional derivative, is a generalization of the classical Fickian diffusion 
equation and describes the asymptotic behavior of a plume of subdiffusive particles in the absence of advection. A 
more general fractional differential equation was derived from CTRW in Refs. [3~l| and [3^]. The power-law kernel 
in 9q7° decays slowly and macroscopically represents the long-range correlations induced by the obstacles that the 
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FIG. 1: Comparison between the analytical result (Eq. (J5J, solid line) and the Monte Carlo variance (dots) as a function of 
time for a = 0.5. The agreement is remarkable. 



particles encounter [32l. |33|. If a > 1, the correlations decay sufficiently quickly so that normal diffusion is recovered 
0- 

In the context of contaminant transport it is of utmost importance to determine the evolution of the plume spread 
as a function of time. This physical quantity corresponds to the variance of the propagator. From the fractional 
diffusion equation it may be expressed in closed form as: 

f) 2 2D 
Because a < 1, the plume migration is subdiffusive, as expected. 



III. HIGHER-ORDER CORRECTIONS TO FRACTIONAL DIFFUSION 



In principle, Eq. ([5|) provides a reliable estimate of the actual contaminant spread when the diffusion limit is 
attained, that is, k « 1 and ut <C 1, so that the fractional diffusion equation is a good asymptotic expansion of 
CTRW. To assess this statement, we compare the analytical variance ([5]) with Monte Carlo simulation results for two 
values of a. Monte Carlo estimates, which are obtained by simulating the microscopic dynamics of the particles, will 
be assumed as a reference solution of CTRW. In the first example, 10 5 particles were followed up to a time £ max = 1, 
with the simulation parameters a — 0.5, a = 0.5, and t = 10~ 4 (see Fig. [1]). In this case the fractional diffusion 
prediction (JSJ) agrees perfectly with the Monte Carlo simulation. 

Next we let a = 0.8, with a = 0.5, r = 10 -4 , and t max = 1 (see Fig. [2]). In this case the fractional diffusion 
prediction underestimates the actual Monte Carlo spread, so that fractional diffusion cannot be considered a good 
approximation of CTRW. The error introduced by the fractional diffusion prediction is far beyond the fluctuations due 
to the finite Monte Carlo statistics. We remark that in both cases the explored time scales are such that r <C t max , 
thus ensuring that the diffusion limit has been attained. 

What is the origin of the discrepancies shown in Fig. [5]? Straightforward but tedious calculations show that the 
expansion of w(u) to second order in u yields: 

w(u) = 1 — Cq,(ut) q — c\{ut) + o(u 2 ), (6) 

where c\ = uj(a — 1). Therefore, when the exponent a is small, the linear term in u in Eq. ^ is expected to play no 
significant role in the limit ut -C 1, thus justifying the accuracy of the fractional diffusion equation predictions. In 
contrast, if a — * 1, the two terms in the expansion ([6]) become comparable and the effects of the linear contribution 
will no longer be negligible (l7ll25j. More precisely, when a < 1 the contribution of c\{ut) will always be subdominant 
with respect to c a (ur) a (for u — > 0), so that after a sufficiently long time it ^> t) we always recover the fractional 
diffusion equation as the correct asymptotic limit of CTRW. In this respect the term C\{ut) introduces pre-asymptotic 
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FIG. 2: Comparison between the analytical (Eq. (|5). solid line) and Monte Carlo results (dots) as a function of time for a — 0.8. 
There is a clear disagreement. 



corrections to the behavior of the subdiffusive particles. However, the values of a and r are usually imposed by 
experiment and inherently linked to the choice of w(t), that is, to the microscopic dynamics of the particles. It 
is therefore interesting to systematically explore the behavior of the CTRW and fractional diffusion solutions for 
intermediate time scales where the diffusion limit ut <C 1 holds but the truncation of (J6j) to the first non-constant 
term might not be appropriate, depending on the value of a. 

To quantitatively assess the relevance of the contribution ci(ur) in the expansion of w(u), we substitute Eq. ([5]) 
into the propagator (Q]), to obtain: 

uP(k, u) - 1 = P(k, u). (7) 

C a T a U a 1 + C\T 

By making use of the formal properties of the Riemann-Liouville operators, Eq. {J} may be finally recast into a 
modified fractional derivative equation 

^P(M) =D a d 1 ^^P(x,t)+q a , 1 d 1 - t a ^ t P(x,t)Y (8) 

where q a ^\ — c\t /c a r a [44j. If a is sufficiently small, we can expand 
1 u 1 -" ( 1 



C\T C a T a \ 1 + q a lit 



,1-a 



i- fc y-« -^-+0^), (9) 



so that we recover the standard fractional diffusion equation. If instead a is of the order of unity, the two terms 
c Q (wr) Q and ci(ur) (which are mirrored in the Z? Q <9g7 a and <?a,i5o7 Q ^ operators, respectively) are in competition, 
and their specific contributions to the overall functional form of the propagator P(x, t) cannot be neglected. We argue 
that Eq. (JSJ) , which is the central result of our work, also provides a suitable description of the particle dynamics for 
pre- asymptotic time scales. To support this argument we compute the spread associated to Eq. ([5]). By definition, 

(x 2 (t)) =C- 1 \% \ 1. (10) 

The inverse C transform appearing in Eq. (|10p can be computed numerically [34], [35[ . The curves thus obtained are 
again compared with Monte Carlo simulation. Figure [3] clearly shows that the modified Eq. ([5]) provides a reliable 
estimate of the system evolution even at pre-asymptotic time scales. 



IV. CONCLUSIONS 



The fractional diffusion equation arises as an asymptotic approximation to an exact master equation and suffers 
from some limitations [20]. As an example, although anomalous diffusion is often experimentally found to be a 
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FIG. 3: Comparison between the variance (|10p as obtained by numerical inversion (solid line) and Monte Carlo simulations 
(dots) for a = 0.8. The two curves are in excellent agreement. 



transient phase, which after a suitable time scale generally relaxes to Fickian diffusion, the fractional diffusion equation 
cannot take into account for this transition because the anomalous behavior is assumed to hold even at infinite time 
@;H, SEl, [HI- Moreover, the CTRW formulation in terms of a generalized Fokker-Planck equation with a memory 
kernel recently proposed e.g. in Q can take into account more general expressions for w(u) (as compared to to those 
allowed by the fractional derivative approach), as well as boundary and initial conditions in multiple dimensions. In 
this respect, the CTRW formalism, thanks to its greater flexibility, can capture and quantify more general tranport 
instances, as shown e.g. in [3rj 

It must be also emphasized that not all anomalous diffusion phenomena can be represented in a fractional derivative 
formalism, since anomalous diffusion is usually context-dependent many specific realizations do not fall within this 
framework (even though enlarged so to include pre-asymptotic corrections). This is the case, for instance, of the 
diffusion problems studied by [38| . where slower than Fickian, but faster than algebraic decays are found, or by [39| . 
where anomalous transport in a fracture network gives rise to a scaling which is at variance with the power-law 
predicted by the fractional derivative framework. 

Nonetheless, a relevant feature of the fractional diffusion equation (in spite of its limitations) is that the fractional 
derivative formulation may easily include external fields in a simple manner and is naturally suitable for solving 
boundary value problems 0, |4l[. In this respect, the fractional diffusion equation has been recently shown to act as 
a unifying framework for the quantitative description of different physical phenomena where anomalous diffusion plays 
a significant role 0, Moreover, many standard mathematical techniques for solving partial differential equations 
are readily applicable to the fractional diffusion equation. 

In the context of particle transport in heterogeneous materials the fractional diffusion formulation may be used 
to analyze the evolution of the contaminant plume spread in time. Comparing the model prediction with measured 
experimental data could then lead to an estimate of the exponent a characterizing the system. It is therefore important 
to quantify the limit of validity of the fractional diffusion asymptotic subset. In this paper we have shown that, if 
a — * 1, the pre-asymptotic corrections to fractional diffusion might significantly affect the model predictions. Their 
neglect would induce gross errors in the model which would distort the estimate of a. To overcome this problem, we 
have derived a modified transport equation involving fractional derivatives: we have shown that the particle spread 
predicted by this model is in excellent agreement with Monte Carlo simulation also at the pre-asymptotic time scales. 
This modified equation might be suitable to explore the subdiffusive dynamics of physical systems close to the limit 
a-»l. 

Finally, we remark that analogous pre-asymptotic corrections would come into play also in the case of diffusive 
dynamics (a > 1), if we were to adopt a power-law decaying distribution for the waiting times. These corrections 
would have a significant role close to the limit a — > 1 and their origin is in the nature of the functional form of w(t). 
Only, the role of the two contributions would be the opposite with respect to the case examined in this work: the 
dominant contribution would come from the term in u, the term in u a being subdominant for u — > 0, as expected. 
As a particular case, it is also possible to choose the (arbitrary) waiting time distribution so that the coefficient C\ 
is identically zero: a widely adopted choice is assuming a Levy stable law for w(t), so that its Laplace transform 
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expansion would not contain the term in u. 
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